{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 作业2：PyTorch基础+线性模型"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 第1题"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(a) 生成一个大小为 $[n \\times p] = [200 \\times 50]$ 的数据矩阵 `x`，用正态分布 N(0, 2) 填充。随机数种子设为123456。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import torch\n",
    "\n",
    "np.random.seed(123456)\n",
    "torch.manual_seed(123456)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "import torch.distributions as D\n",
    "import math\n",
    "n = 200\n",
    "p = 50\n",
    "# 完成后续的代码\n",
    "norm1 = D.Normal(loc=torch.tensor([0.0]), scale=torch.tensor([math.sqrt(2.0)]))\n",
    "x = norm1.sample((n,p))\n",
    "x = torch.squeeze(x)\n",
    "\n",
    "# 检查 x 的大小，方便 debug\n",
    "assert x.shape == (n, p), \"x 形状有误\"\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(b) 生成一个长度为 $p$ 的向量 `beta`，每个元素服从均匀分布 Uniform(-1, 1)。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "uniform_dist = D.Uniform(-1, 1)\n",
    "beta = uniform_dist.sample((p,))\n",
    "\n",
    "# 检查 beta 的长度，方便 debug\n",
    "assert beta.shape == (p,), \"beta 长度有误\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(c) 生成一个长度为 $n$ 的向量 `eps`，每个元素服从独立正态 $N(0, 0.1)$。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "norm2=D.Normal(loc=torch.tensor([0.0]),scale=torch.tensor([math.sqrt(0.1)]))\n",
    "eps = norm2.sample((n,))\n",
    "eps = torch.squeeze(eps)\n",
    "# 检查 eps 的长度，方便 debug\n",
    "assert eps.shape == (n,), \"eps 长度有误\"\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(d) 创建向量 `y`，令其在数学上等于 $y=X\\beta+\\epsilon$。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "y = torch.matmul(x,beta)+eps\n",
    "\n",
    "# 检查 y 的长度，方便 debug\n",
    "assert y.shape == (n,), \"y 长度有误\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(e) 回归问题：给定数据 `x` 和 `y`，估计 `beta` 的取值。以 MSE 为损失函数，编写 Python 函数 `loss_fn_reg(bhat, x, y)`，用来返回任意 $\\hat{\\beta}$ 下的目标函数值。请用基础的矩阵和向量运算实现。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "def loss_fn_reg(bhat, x, y):\n",
    "    yhat = torch.matmul(x, bhat)\n",
    "    mse = torch.mean((y-yhat)**2)\n",
    "    return mse\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(f) Pytorch 中也提供了 MSE 损失函数，参见[其文档](https://pytorch.org/docs/stable/generated/torch.nn.MSELoss.html#torch.nn.MSELoss)。其用法是先建立一个损失函数对象，然后将 $\\hat{y}$ 和 $y$ 作为参数传入。请利用这种方法计算如下给定 $\\hat{\\beta}$ 后的损失函数值，并与你自己的函数结果进行对比。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "tensor(145.7161)\n",
      "tensor(145.7161)\n"
     ]
    }
   ],
   "source": [
    "import torch.nn as nn\n",
    "\n",
    "bhat = torch.ones(p)\n",
    "yhat = torch.matmul(x,bhat)  # 替换这里的代码计算 yhat\n",
    "\n",
    "mse_reg = nn.MSELoss()\n",
    "loss1 = mse_reg(yhat, y)\n",
    "\n",
    "loss2 = loss_fn_reg(bhat, x, y)\n",
    "\n",
    "print(loss1)\n",
    "print(loss2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "【说明文字】两种方式计算的损失值应该是等价的，因为它们都用于衡量模型的预测与实际标签之间的均方误差，所以两种方法的函数结果一致"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(g) 利用 PyTorch 的自动微分功能，计算 MSE 损失函数在 $\\hat{\\beta}=\\mathbf{1}_p$ 处的梯度，其中 $\\mathbf{1}_p$ 是一个元素全为1的向量。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 在此处插入代码"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 第2题"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(a) 与第1题类似创建变量 `x` 和 `beta`，但使用不同的 `n` 和 `p`。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(123456)\n",
    "torch.manual_seed(123456)\n",
    "n = 150\n",
    "p = 6\n",
    "x = norm1.sample((n,p))\n",
    "x = torch.squeeze(x)\n",
    "beta = uniform_dist.sample((p,))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(b) 定义函数 `sigmoid(x)`，其中 `x` 是一个 Tensor，$\\mathrm{sigmoid}(x)=e^x/(1+e^x)$。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sigmoid(x):\n",
    "    \n",
    "    return torch.exp(x) / (1+torch.exp(x))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(c) 根据分布 $Y|X\\sim Bernoulli(\\mathrm{sigmoid}(X\\beta))$，生成 $Y$ 的随机数。提示：参照1.4节的方法，先计算 Bernoulli 分布的参数**向量**，然后生成随机数。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "tensor([1., 0., 1., 0., 0., 1., 1., 0., 0., 1., 1., 0., 1., 0., 0., 0., 0., 0.,\n",
      "        1., 0., 0., 1., 0., 0., 1., 1., 1., 0., 0., 0., 0., 0., 0., 1., 0., 0.,\n",
      "        0., 0., 0., 1., 0., 1., 1., 1., 1., 0., 1., 0., 0., 0., 1., 0., 1., 0.,\n",
      "        0., 0., 1., 0., 0., 0., 1., 1., 0., 0., 0., 1., 1., 0., 0., 0., 1., 1.,\n",
      "        0., 1., 0., 0., 1., 1., 1., 1., 1., 0., 0., 0., 0., 0., 1., 1., 1., 0.,\n",
      "        0., 1., 1., 1., 1., 1., 0., 0., 0., 1., 0., 1., 0., 1., 1., 1., 0., 1.,\n",
      "        0., 0., 1., 0., 1., 0., 0., 1., 1., 1., 0., 0., 0., 1., 1., 1., 0., 1.,\n",
      "        0., 1., 0., 0., 1., 0., 0., 1., 1., 1., 0., 0., 0., 0., 0., 0., 1., 1.,\n",
      "        1., 1., 0., 0., 0., 0.])\n"
     ]
    }
   ],
   "source": [
    "p_= sigmoid(torch.matmul(x,beta))\n",
    "bern = D.Bernoulli(p_)\n",
    "y = bern.sample()\n",
    "print(y)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(d) 已知 $Bernoulli(\\rho)$ 分布的对数密度函数为 $\\log p(y;\\rho)=y\\log \\rho + (1-y) \\cdot \\log(1-\\rho)$。根据此信息，推导出给定 $\\hat{\\beta}$ 时的对数似然函数，并编写损失函数 `loss_fn_logistic(bhat, x, y)`，返回**负**对数似然值。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "【说明文字】"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [],
   "source": [
    "def loss_fn_logistic(bhat, x, y):\n",
    "    phat = sigmoid(torch.matmul(x,bhat))\n",
    "    #计算对数似然函数\n",
    "    log_like = torch.mean(y * torch.log(phat) + (1 - y) * torch.log(1 - phat))\n",
    "    ne = -log_like\n",
    "    return ne"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(e) Pytorch 中也提供了 BCELoss 损失函数，参见[其文档](https://pytorch.org/docs/stable/generated/torch.nn.BCELoss.html)。其用法是先建立一个损失函数对象，然后将 $\\hat{\\rho}$ 和 $y$ 作为参数传入。请利用这种方法计算如下给定 $\\hat{\\beta}$ 后的损失函数值，并与你自己的函数结果进行对比。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "tensor(1.3206)\n",
      "tensor(1.3206)\n"
     ]
    }
   ],
   "source": [
    "bhat = torch.ones(p)\n",
    "rhohat = sigmoid(torch.matmul(x,bhat)) # 替换这里的代码计算 rhohat = sigmoid(x * betahat)\n",
    "\n",
    "bce_logistic = nn.BCELoss()\n",
    "loss1 = bce_logistic(rhohat, y)\n",
    "\n",
    "loss2 = loss_fn_logistic(bhat, x, y)\n",
    "\n",
    "print(loss1)\n",
    "print(loss2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "【说明文字】】 结果一样的原因是 bce_logistic 和 loss_fn_logistic 这两个函数计算的损失值实际上是相同的。这是因为它们两个分别实现了相同的二分类的二元交叉熵损失函数，只是使用了不同的计算方式。\n",
    "\n",
    "bce_logistic 使用了 PyTorch 提供的 nn.BCELoss() 函数，这是一个高度优化的损失函数，专门用于二分类问题。它接受模型的预测概率 rhohat 和实际标签 y，并计算二元交叉熵损失。\n",
    "\n",
    "loss_fn_logistic 函数是自定义的损失函数，它也计算了二元交叉熵损失，但是是基于模型的预测概率 phat 和实际标签 y 计算的，使用了手动编写的计算方式。\n",
    "\n",
    "由于这两个函数都实现了相同的损失计算逻辑，因此它们的结果是一样的。这是因为这个问题是一个简单的二分类问题，所以无论使用哪种方式计算损失，都应该得到相同的结果。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(f) 利用 PyTorch 的自动微分功能，计算上述损失函数在 $\\hat{\\beta}=\\mathbf{0}_p$ 处的梯度，其中 $\\mathbf{0}_p$ 是一个元素全为0的向量。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 在此处插入代码"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 第3题"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(a) 多分类问题的数据通常包括数据阵 $X$ 和标签向量 $l$，其中标签为整数。在计算损失函数时，我们需要先将 $l$ 转换成多项分布的0-1数据，即所谓 One-hot 编码。运行并观察下面的代码。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "tensor([1, 2, 2, 1, 0, 3, 3, 3, 3, 0, 3, 0, 0, 2, 2, 0, 3, 0, 3, 3])\n",
      "torch.Size([200, 4])\n",
      "tensor([[0, 1, 0, 0],\n",
      "        [0, 0, 1, 0],\n",
      "        [0, 0, 1, 0],\n",
      "        [0, 1, 0, 0],\n",
      "        [1, 0, 0, 0],\n",
      "        [0, 0, 0, 1],\n",
      "        [0, 0, 0, 1],\n",
      "        [0, 0, 0, 1],\n",
      "        [0, 0, 0, 1],\n",
      "        [1, 0, 0, 0]])\n"
     ]
    }
   ],
   "source": [
    "np.random.seed(123456)\n",
    "torch.manual_seed(123456)\n",
    "n = 200  # 样本量\n",
    "p = 10   # 变量数\n",
    "k = 4    # 类别数\n",
    "x = torch.randn(n, p)\n",
    "l = torch.tensor(np.random.choice(range(4), size=n, replace=True), dtype=int)\n",
    "print(l[:20])\n",
    "\n",
    "y = torch.nn.functional.one_hot(l)\n",
    "print(y.shape)\n",
    "print(y[:10])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(b) 创建矩阵 `W`，大小为 $k \\times p$，用 N(0, 1) 填充其取值。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "w = torch.tensor(np.random.normal(loc=0,scale=1,size=(k,p)))\n",
    "\n",
    "# 检查 w 的形状，方便 debug\n",
    "assert w.shape == (k, p), \"w 形状有误\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(c) 接下来计算对 $Y$ 的概率预测值，其中每个 $Y_i$ 观测对应一个等长的概率向量 $p_i$，而 $p_i=\\mathrm{Softmax}(Wx_i)$。首先计算 $Wx_i$，其中 $x_i$ 是第 $i$ 个观测。由于 $X$ 是把 $x_i$ 按行组合，因此矩阵形式表达为 $U=XW'$，其中 $U$ 的第 $i$ 行即为 $Wx_i$。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "u = torch.matmul(x.float(), w.t().float())  # 替换这里的代码\n",
    "\n",
    "# 检查 u 的形状，方便 debug\n",
    "assert u.shape == (n, k), \"u 形状有误\"\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们先测试一下 $\\mathrm{Softmax}(Wx_{100})$ 的结果，观察其元素和是否为1。代码中的 `dim=0` 意思是对第一个下标方向计算 Softmax，由于 `U[99]` 是一个向量，因此第一个下标方向就是该向量本身。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "tensor([1.0882e-05, 2.1499e-04, 9.9927e-01, 5.0517e-04])"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "torch.softmax(u[99], dim=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "而为了对 $U$ 的每一行都计算 Softmax，我们可以直接对整个 `U` 矩阵用 `torch.softmax`，其中 `dim` 需指定为1，意思是对第二个下标方向求 Softmax，即对 $U$ 的每一行。原理类似于1.3节的按坐标轴汇总。请完成该计算，得到矩阵 $P$，其中 $P$ 的第 $i$ 行即为 $p_i$。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "p = torch.softmax(u, dim=1)  # 替换这里的代码\n",
    "\n",
    "# 检查 p 的形状，方便 debug\n",
    "assert p.shape == (n, k), \"p 形状有误\"\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(d) 根据 `y` 和 `P` 两个矩阵，即可根据公式得到对数似然函数值。总结上述步骤，编写损失函数 `loss_fn_softmax(w, x, y)`，返回**负**对数似然值。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "def loss_fn_softmax(w, x, y):\n",
    "    # 编写函数主体，替换这里的代码\n",
    "    # 计算预测概率分布\n",
    "    p = torch.softmax(u, dim=1)\n",
    "\n",
    "    # 计算负对数似然值\n",
    "    loss = -torch.sum(y * torch.log(p))/x.shape[0]\n",
    "\n",
    "    return loss"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(e) Pytorch 中也提供了 CrossEntropyLoss 损失函数，参见[其文档](https://pytorch.org/docs/stable/generated/torch.nn.CrossEntropyLoss.html)。其用法是先建立一个损失函数对象，然后将 $U$ 和 $l$ 作为参数传入（注意 $U$ 是经过 Softmax **之前**的矩阵，$l$ 是**原始**的标签）。请利用这种方法计算损失函数值，并与你自己的函数结果进行对比。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "tensor(3.9089)\n",
      "tensor(3.9089)\n"
     ]
    }
   ],
   "source": [
    "ce_softmax = nn.CrossEntropyLoss()\n",
    "loss1 = ce_softmax(u, l)\n",
    "\n",
    "loss2 = loss_fn_softmax(w, x, y)\n",
    "\n",
    "print(loss1)\n",
    "print(loss2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "【说明文字】结果一样的原因是 nn.CrossEntropyLoss() 和 loss_fn_softmax 这两个函数计算的损失值实际上是相同的，但它们的输入和计算方式略有不同。\n",
    "\n",
    "nn.CrossEntropyLoss() 函数是 PyTorch 提供的用于多分类问题的损失函数。它要求标签 l 是包含每个样本对应类别的索引的向量，而不是 One-hot 编码的形式。它使用未经 Softmax 处理的模型输出 u，并将其与类别索引 l 进行比较，计算交叉熵损失。\n",
    "\n",
    "loss_fn_softmax 函数是自定义的损失函数，它用于计算多分类问题的负对数似然损失。它使用了 Softmax 处理后的模型输出 w，并与 One-hot 编码的标签 y 进行比较，计算负对数似然损失，最后取平均损失。\n",
    "\n",
    "虽然这两个函数的输入和计算方式不同，但由于它们都是用于多分类问题的标准损失函数，所以在这个问题中，它们计算的损失值应该是相同的。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(f) 利用 PyTorch 的自动微分功能，计算上述损失函数在 $W=O$ 处的梯度，其中 $O$ 是一个元素全为0的矩阵。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 在此处插入代码"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
